Magnetic and Transport Properties of a Coupled Hubbard Bilayer with Electron and 

Hole Doping 
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The single band, two dimensional Hubbard Hamiltonian has been extensively studied as a model 
for high temperature superconductivity. While Quantum Monte Carlo simulations within the dy- 
namic cluster approximation are now providing considerable evidence for a d-wave superconducting 
state at low temperature, such a transition remains well out of reach of finite lattice simulations 
because of the "sign problem". We show here that a bilayer Hubbard model, in which one layer 
is electron doped and one layer is hole doped, can be studied to lower temperatures and exhibits 
an interesting signal of d-wave pairing. The results of our simulations bear resemblance to a recent 
report on the magnetic and superconducting properties of Ba2Ca3Cu40sF2 which contains both 
electron and hole doped Cu02 planes. We also explore the phase diagram of bilayer models in 
which each sheet is at half-filling. 

PACS numbers: 



I. INTRODUCTION 

The single band, two dimensional Hubbard Hamilto- 
nian provides one possible microscopic model for pairing 
which is driven by electronic correlations rather than the 
interactions of electrons with the lattice. Many analytic 
and numeric^ treatments suggest that there may indeed 
be a superconducting phase at low temperature away 
from half-filling in this model. The issue is a difficult 
one, however, owing to the likely existence of a variety 
of different phases which are close in energy on the one 
hand, and the nature of the approximations made in the 
solution on the other. Exact diagonalization studies^, 
while very useful, are typically on lattices of only a few 
tens of sites, and hence finite size effects are a consider- 
able concern. Quantum Monte Carlo (QMC^, which 
can in principle address the issue in an unbiased way (on 
lattices an order of magnitude or more larger than di- 
agonalization) has been unable to access sufficiently low 
temperatures due to the 'minus sign problem—. 

Recently, progress has been made using improved nu- 
merical methods. The 'density matrix renormalization 
group' has pushed forward from one dimension to ad- 
dress geometries of many coupled chains^. The dy- 
namic cluster approximation has improved on dynami- 
cal mean field treatments by showing the robustness of 
a finite temperature transition to a superconducing state 
as an increasingly fine momentum grid is incorporated 
in the self-energy^. Nevertheless, there is still numeric 
work which contests the conclusion that the two dimen- 
sional Hubbard Hamiltonian has long range d-wave pair 
correlations^. 

In this paper, we present determinant Quantum Monte 
Carlo (DQMC) calculations of a bilayer Hubbard model 
for which we are able to attain much lower temperatures 
than the single layer case. Specifically, by doping the 
two layers symmetrically about half-filling, p = 1, we 
find that the sign problem is greatly reduced, allowing 
simulations at temperatures which are roughly two or- 



ders of magnitude below the bandwidth, T w VF/100. 
In single layer simulations of the doped system, the low- 
est attainable temperatures are T « IF/40. Previous 
DQMC studies of bilayer models have looked at the case 
when both layers are half-filled, and examined magnetic 
order-disorder transitions which occur as the interlaycr 
hopping is increased^. A decreasing interlayer hopping 
monotonically reduces the pairing correlations in this sit- 
uation. 

Our work is partially motivated by studies of 
Ba2Ca3Cu40gF2 and Ba2Ca2Cu40eF2 which are an ex- 
perimental realization of materials in which electron and 
hole doped sheets coexist within the family of cuprate 
superconductors^. In the former, four-layered com- 
pound, the two outer planes are electron-doped with 
N e 0.06 — 0.08, while the two inner planes are hole 
doped roughly symetrically, that is Nh ~ 0.06 — 0.08. 
The superconducting transition temperature is T c = 55 
K, and pairing coexists with long range antiferromag- 
netic order with Neel temperature TV = 100 K. The 
latter, three-layered compound has outer plane doping 
N e W 0.06 — 0.08, but a larger inner plane doping 
Nh ~ 0.13. Its superconducting T c = 76K with only 
short range antiferromagnetic correlations. This is at- 
tributed to a decoupling of the magnetism of the elec- 
tron doped outer planes by the large doping of the inner 
planoii. 



II. MODEL AND METHODOLOGY 

In order to model such materials, we consider the two 
layer Hubbard Hamiltonian, 

H = -t ( c j» c i» + k ') 
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The first term is the usual hopping of electrons between 
near neighbor sites i and j of a two dimensional square 
lattice. Unless otherwise stated, the results in this pa- 
per are for two coupled 8x8 lattices. The electrons in 
the kinetic energy term have a spin index a =T, j and 
also a layer index m = 1,2. The second term is an in- 
tcrlaycr hopping. The third term is a layer-dependent 
chemical potential. We will choose fii = —/j,2 to produce 
layers which have opposite dopings. Finally, electrons of 
opposite spin on the same site of the same layer feel a 
repulsion U. 

Our simulations employ the DQMC algorithm^^ in 
which a path integral is written for the partition func- 
tion, the fermion interactions are replaced by a coupling 
to an auxiliary Hubbard-Stratonovich field, and then the 
fermion degrees of freedom are integrated out analyti- 
cally. The method produces exact results on the lattice 
sizes considered, apart from 'Trotter' errors associated 
with the imaginary time discretization, which we have 
verified are smaller than our statistical error bars. 

The magnetic properties are determined from the spin- 
spin correlations, 



,( jmt n jml , (2) 

which arc independent of layer index m because of our 
choice of symmetric doping and the particle-hole symme- 
try of the Hubbard Hamiltonian. The Fourier transform 
gives the structure factor, 



(3) 



At half-filling, S(q) is largest at the antiferromagnctic 
wavevector q = (tt, 7t, it). 

A first insight into the metal-insulator transition can 
be obtained from the zero momentum spectral function 
(density of states) A(u>) which is determined from the 
Greens function, 

Gi_j, mCT (r) = (c iroo -(r)c] roCT (0)) 
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using the maximum entropy method^. 

The dc conductivity a^ c also characterizes the metal- 
insulator transition, and is measured from the current- 
current correlation function, 

j x (l,r) = e Hr jS,0)e- Hr 

jxO-,0) = ^yi^H-^rnjClmg ~ c l ma c \+xma) 
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This imaginary time quantity, which comes directly out 
of the determinant QMC simulations, is related to the 
real frequency response by the fluctuation-dissipation 
theorem, 



A ra (q; r) 
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As discussed in [14j, at sufficiently low temperatures we 
can replace ImA by its low frequency behavior ImA « 
CTdcW, leading to the relation, 
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This enables us to obtain the conductivity directly from 
the imaginary time response without the necessity for an- 
alytic continuation, which is more difficult for two par- 
ticle response functions, like the current-current correla- 
tor, than for the single particle Greens function, owing 
to their larger fluctuations. 

To describe superconductivity, we compute the corre- 
lated pair field susceptibility, P a , in different symmetry 
channels, 
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The correlated susceptibility P a takes the expectation 
value of the product of the four fermion operators enter- 
ing Eq. [51 We also define the uncorrelatcd pair field sus- 
ceptibility P a which instead computes the expectation 
values of pairs of operators prior to taking the product. 
Thus, for example, in the s-wave channel, 

Ps = rfr(c u (r) CiT (r) ^(0)^(0)) 
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P * = / dT(c n (T)cl(0)) <q T (r)c] T (0)> (9) 

P a includes both the renormalization of the propagation 
of the individual fcrmions as well as the interaction ver- 
tex between them, whereas P a includes only the former 
effect. Indeed by evaluating both P and P we are able 
to extract^ the interaction vertex T, 
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If T a P a < 0, the associated pairing interaction is attrac- 
tive. In fact, rewriting Eq. 10 as, 



P a = 



i + r a P a 
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FIG. 1: (Color online) Phase diagram for the half-filled bilayer 
Hubbard model. A paramagnetic metallic phase is present at 
weak coupling. At large coupling there is a transition from an 
antiferromagnetic Mott-insulating phase to a paramagnetic 
band-insulating phase. The phase boundaries obtained by 
the conductivity a and density of states at the Fermi level, 
^4(0), are consistent. 

suggests that T a P a — > —1 signals a superconducting in- 
stability. We will discuss this criterion in more detail in 
the coming sections. 

III. BILAYER PHASE DIAGRAM AT 
HALF-FILLING 

We begin with the phase diagram at half-filling, that is 
when /ii = fi2 = 0, and both layers have equal occupation 
pi = p2 = 1. (Note there is no sign problem in this case 
because of particle- hole symmetry.) Here we do not ex- 
pect superconductivity. Nevertheless there is an interest- 
ing competition between Mott insulating behavior when 
U is the dominant energy scale, and band insulating be- 
havior for large tj_. Indeed, increased interlayer coupling 
suppresses the antiferromagnetic correlations which are 
present in the Mott phase, since t± promotes the forma- 
tion of interlayer singlets between the two spatial sites 
immediately above and below each other. These spin- 
singlets arc magnetically decoupled, destroying long 
range spin order. Earlier determinant QMC studies de- 
termined the critical value of t± 1.6 for this AF-PM 
transition^. 

The strong coupling region of Fig. [T] exhibits this phe- 
nomena, and yields a (tj_/t) c consistent with the ear- 
lier study-i^. At weak coupling, however, this insulator- 
insulator transition is replaced by a metallic phase. Pre- 
vious cluster DMFT— studies of the bilayer model show 
a phase diagram which is in qualitative agreement with 
Fig. [TJ We will compare the results of the two methods 
in more detail at the end of this section. First, we will 
describe in detail how this phase diagram is obtained. 

In Fig. [5] the density of states at the Fermi surface, 
A(u> = 0) is shown for four temperatures along a hori- 
zontal cut through the phase diagram at fixed t±/t = 2. 
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FIG. 2: (Color online) Density of states at the Fermi surface 
A(0) for t± = xx. At weak U, A(0) rises as T = 1/(3 is 
lowered, indicating a metallic phase with nonzero Fermi level 
density of states. In contrast, at large U, A(0) falls with 
decreasing T, indicating insulating behavior. (U/t) c ~ 2.8. 

At weak coupling, the low temperature limit is non-zero, 
indicating a metallic phase, while at strong coupling, 
A(u> = 0) decreases as T is lowered. We conclude that at 
the crossing point U/t ~ 2.8 a metal-insulator transition 
occurs. 

In Fig. [3] we see that the conductivity a^c similarly 
can determine the location of the metal-insulator phase 
boundary. Here a change in the temperature behavior of 
the conductivity, from increasing as T is lowered (metal- 
lic) to decreasing when T is lowered (insulating) occurs 
at U/t « 2.6 when the interlayer hopping is t±/t = 3.4. 

Multiple horizontal (constant t±/t) cuts through the 
phase diagram similar to those of Figs. [2H2] were used to 
generate the metal-insulator phase boundary of Fig. [JJ 
Note the consistency of the locations of the critical inter- 
action strengths between those obtained from the density 
of states A(uj = 0) (red squares in Fig. [JJ and the con- 
ductivity Ode (green diamonds in Fig.[TJ). 

In this bilayer model, at half- filling fi± = /12 = 0, the 
suppression of the zero frequency spectral weight can 
come from any of three mechanisms: the opening of a 
band gap at sufficiently large t±, a "Slater gap" created 
by antiferromagnetic fluctuations which can form on a 
scale set by the exchange constant J oc t 2 /U, and a "Mott 
gap" between the upper and lower Hubbard bands when 
U exceeds the bandwidth W. (The bandwidth W = 8t 
at t± =0.) In general, these different insulating phases 
cross over to each other more or less smoothly, although 
the Slater insulator can be distinguished by the presence 
of long range spin correlations. Fig. [4] shows the full fre- 
quency dependence of the density of states at U/t = 3 
and three values of t±, all of which exhibit a gap in A(ui). 
(The non-zero residual values of A(u>) for t± = 1.4 and 
4.0 will be driven to zero if (3 is increased. See Fig. [2]) 
From Fig. 0] we infer that the phase diagram is insulating 
all along the vertical line U/t = 3 in Fig. [TJ 

In contrast, Fig. [5J which shows the same three values 
of ij_ except at weaker coupling, U/t = 2, clearly ex- 
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FIG. 3: (Color online) The conductivity a^ c for a horizontal 
cut (fixed t±/t = 3.4) and varying U/t through the phase dia- 
gram. Values at four inverse temperatures are given. As with 
the density of states at the Fermi energy, A(u> = 0), shown 
in Fig. [2] the conductivity exhibits a crossing pattern which 
gives the location of the metal-insulator phase boundary: cr<j c 
increases as f3 increases (metallic behavior) below U/t ~ 2.6, 
and falls as j3 increases above this value. 




FIG. 5: (Color online) The same as Fig. |4j except U/t = 2. 
Although t± — 0.0 and 4.0 are still insulating, the density of 
states for t± = 1.4 has a peak at w = and is metallic in 
character. 
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FIG. 6: (Color online) Real space spin correlations at U/t = 5. 
As t± increases, the antiferromagnetism is suppressed. The 
inverse temperature (3 = 14. 



FIG. 4: (Color online) Density of states A(lu) at U/t = 3 
and inverse temperature /3 = 14, showing insulating behavior 
at all values of interlayer coupling. t± = and t± = 1.4 
are Mott/Slater insulating phases with a gap produced by 
a combination of the on-site repulsion and antiferromagnetic 
spin correlations. t± = 4.0 has a gap which is primarily band 
insulator in character. 



hibits metallic behavior for the intermediate value of the 
interlayer hopping. This is one indication of the outward 
extent of the metallic region from U/t = in the phase 
diagram of Fig. [T] 

We turn now to the spin correlations. Fig. [S] shows 
the real space spin correlations for U/t = 5 and different 
interlayer hoppings. t± drives the formation of interlayer 
singlets which interfere with the magnetic order. A finite 
size scaling analysis is shown in Fig. [7] where the struc- 
ture factor is plotted as a function of the inverse linear 
system size. Spin wave theory predicts^ that the finite 
size corrections to S(ir, 7t, it) should be linear in 1/N X , 
with the N x — > oo intercept proportional to the square of 
the order parameter. We see that the order parameter is 



non-zero for tx/t = 0.6 and 1.4 and is zero for tx/t = 2.8 
and 3.4. Somewhere in the vicinity of t±/t f=s 2, the long 
range magnetic order is destroyed. Fig. [5] shows a similar 
finite size scaling analysis for weaker coupling, U/t = 2. 
There is no long range magnetic order for any value of 
tx/t. 

Multiple vertical (constant U /t) cuts through the 
phase diagram similar to those of Figs. were used to 
generate the limit of the antiferromagnetically ordered re- 
gions of the phase diagram Fig.[TJ This value is consistent 
with previous DQMC studies^ and cluster DMFTi£. 

We conclude this section with a more quantitative com- 
parison of Fig. [T] with the results obtained in cluster 
DMFT— . At strong coupling, the AF insulator to para- 
magnetic (bond) insulator transition is found by both 
methods to have the same value tx/t = 2. Likewise, in 
both approaches, the base of the metallic phase at U = 
extends from tx/t = to tx/t = 4, as indeed it must an- 
alytically from the non-interacting dispersion which has 
bonding and anti-bonding bands, 

ei(k) = — tx + 2i ( cos k x + cos k y ) 
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FIG. 7: (Color online) Scaling of the antiferromagnetic struc- 
ture factor at U — xx. If there are long range correlations, 
5(7r,7r,7r) should grow linearly with lattice size N, so that 
S(n,n , 7r) /N approaches a constant for large N. Spin wave 
theory predicts a 1/N X correction, where N x is the linear lat- 
tice size (N% = N) . Here we see long range order for the three 
smallest values t±/t = 0.4, 1.4, 2.0, but not for the two largest 
values t±/t = 2.8,3.4. 
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FIG. 8: (Color online) Same as Fig. [7] except U/t = 2. 



e 2 (k) = +t ± + 2t ( cos k x + cos k y ) 



(12) 



The extent of the metallic phase as U increases from 
the non-interacting limit differs quantitatively in the two 
methods. The DQMC results reported here indicate an 
upper limit of U/t w 3, while within cluster DMFT the 
metallic region extends out to U/t « 8. The precise ori- 
gin of this disagreement is not clear. The peak of the clus- 
ter DMFT metallic lobe follows the emerging AF-band 
insulator line rather narrowly, and it is possible DQMC 
cannot resolve this small region adequately. While the 
results of Figs. [2] and |4] seem unambiguously to rule out 
metallic behavior much beyond U/t ~ 3, they are on lat- 
tices of finite extent (A=8x8). Cluster DMFT works in 
the thermodynamic limit and hence typically produces 
sharper transitions which can distinguish narrow regions 
of phase space. On the other hand, DQMC incorporates 
the full momentum dependence of the self-energy, in con- 



trast to the 2x2 momentum grid used in cluster DMFT. 



IV. SUPERCONDUCTIVITY IN THE DOPED 
SYSTEM 

Fig.[5]shows a central result of our paper. The product 
of the d-wave pairing vertex and the uncorrelated suscep- 
tibility, TdPd, is seen to turn sharply negative (attractive) 
as the temperature T is lowered. As described in Eq. 11, 
TdPd — 1 in principle would signal a superconducting 
instability. For p w 0.87, TdPd ~ —0.9 . In comparison, 
the most negative TdPd reported^ for the single band 
model is T d P d = -0.45 at half-filling, p = 1.000, and 
T d Pd = -0.25 for doping to p = 0.875. It should be kept 
in mind, however, that the lowest accessible temperature 
in the latter case is f3 = 6/t. At the same (3 = 6/t and 
doping p = 0.875, as seen in Fig. [SI the bilayer system 
has a somewhat more negative TdPd = —0.31. Thus the 
approach of TdPd to —1 seen in the bilayer system is due 
both to a more attractive pairing vertex, but also due to 
the ability of the simulation to reach much colder tem- 
peratures. 

Although wc find the vertex TdPd approaches -1, this 
criterion for an instability is incomplete. One also needs 
to require that the uncorrelated susceptibility P remain 
finite at the transition point. Especially in the situation 
where there is competing order (e.g. antiferromagnetism 
and pairing) it is possible for the uncorrelated suscepti- 
bility of one type of order to be driven to small values 
by the other order, so that even though the vertex ap- 
proaches -1, order in this channel is usurped. Fig. [TOl 
addresses this issue for the bilayer model. Despite the 
fact that TdPd is getting close to —1, the correlated ver- 
tex Pd grows relatively slowly as T is decreased. The 
reason is clear from Fig. [10] in which it is seen that the 
uncorrelated susceptibility is rapidly dropping as T is 
lowered. This is rather different from the doped single 
layer model, where Pd grows as T is lowered. (At half- 
filling in the single layer model Pd declines slightly as T 
is decreased, as found here also in the bilayer model.) 

An interesting feature of Figs. [§] and [TU] is that the 
d-wave attraction is maximal at p w 0.87, whether mea- 
sured via the vertex or the correlated susceptibility. This 
point is made more concretely in Fig. 111! The behavior 
of the d-wave superconducting vertex bears an interest- 
ing resemblance to the superconducting "domes" of the 
cupratc materials in which the transition temperatures 
are maximized a finite distance away from "half-filling" 
(one hole per Cu). Indeed, even the values of the dop- 
ing which maximizes T c and the width of the base of the 
dome are in reasonable quantitative agreement. It is to 
be emphasized that, within the same DQMC method- 
ology, the single layer Hubbard model has a maximum 
pairing vertex at half-filling. Fig. [Tl] also indicates that, 
within the parameter range accessible, the degree of en- 
hancement increases as t± decreases. Eventually we ex- 
pect this trend to reverse, since at t± = 0, the single layer 



6 



.-o 



ICL 

u 



-0.2 
-0.4 
-0.6 
-0.8 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 








•-•p=1.0 




■-.p=0.96 




~p=0.916 
^p=0.869 




«-«p=0.821 













0.05 



0.1 



0.15 



0.2 



0.25 



FIG. 9: (Color online) d-wave pairing vertex as a function 
of temperature for two 8x8 bilayers with interlayer hopping 
t± = 0.6t. The on-site interaction U — 3t. Three fillings are 
shown. Note the close approach to VP = —1, the onset point 
of a pairing instability, and the non-monotonic dependence 
on filling. The greatest tendency to pairing is at p ~ 0.87. 



model, there is a lesser tendency for pairing. (We cannot 
accumulate data for smaller values of t± because the sign 
problem prevents simulations at as low a temperature as 
for the data shown.) 

We turn now to the magnetic properties of the 
doped system and in particular their connection to 
those observed in the cuprate superconductors. Fig. Q2] 
shows the real space spin-spin correlations for p = 
1.00,0.96,0.92,0.87 and 0.82 at = 8/t, U = M, 
T±_ = 0.6. These results have a quantitative similar- 
ity to the Ba2Ca3Cu40gF2 and Ba2Ca2Cu40gF2 mate- 
rials in that the robust magnetic correlations present for 
p = 1.00 and p = 0.96 are dramatically suppressed for 
p = 0.87. A finite Necl temperature T/v is present for the 
four layer compound Ba2Ca2Cu40gF2 which has electron 
and hole dopings N e , Nh ~ 0.06 and absent for the three 
layer compound Ba2Ca2Cu40gF2 which has hole doping 
Nh ~ 0.14 in the central layer. 

Why is the sign problem ameliorated in these bilayer 
simulations? In DQMC for the single layer Hamiltonian, 
the operator couples to the Hubbard Stratonovich 
field hf^- shifted by the chemical potential hi — p. Mean- 
while, riii couples to —hi — p.. At half-filling, fx = 0, 
particle-hole symmetry is reflected in the fact that the 
up and down species couple to the quantities ±/ij which 
are symmetric about zero. The up and down determi- 
nants can be shown to have the same sign, and hence 
their product is positive. For p ^ this symmetry and 
the associated connection between the signs of the two 
determinants is broken, and a sign problem results. (Note 
that for the attractive Hubbard Hamiltonian and 
both couple to hi — p and the two determinants are equal 
at all fillings.) 

Consider now the bilayer system. We have a Hubbard- 
Stratonovich field for each layer. The operators rijif cou- 
ple to hn — p, while mil cou ple to —hn — fx, and n^t cou- 
ple to hi2 + p, and finally n^j, couple to —hi2 + fJ., where 




FIG. 10: (Color online) Correlated (closed symbols) and un- 
corrected (open symbols) d-wave pairing susceptibility as 
a function of temperature for two 8x8 bilayers with inter- 
layer hopping t± — 0.6t. The on-site interaction U = 3t. 
Five fillings are shown. In all cases the vertex is attractive, 
ie. Pd > Pd- The degree of attraction is non-monotonic, first 
increasing with doping, but then declining. 
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FIG. 11: (Color online) d-wave pairing vertex as a func- 
tion of filling for two 8x8 bilayers with interlayer hopping 
t±/t = 0.6,0.7,0.8. The on-site interaction U = 3t and in- 
verse temperature (3 — 14. The greatest tendency to pairing 
is at p w 0.87. 
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FIG. 12: (Color online) Real space spin correlations. At half 
filling, p = 1.00, and for small dopings, p = 0.96, there is a 
strong oscillatory pattern indicative of long range magnetic 
ordepi^. For larger dopings, the spin correlations are sharply 
curtailed. 
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we have explicitly set \i\ = —^2 = — M- What we observe 
is that, to the extent that the Hubbard-Stratonovich vari- 
ables on the two layers are equal, rijiT an d n.j2| are sym- 
metrically coupled about zero. It is possible that this 
tends to lead to a positive determinant product similar 
to the single layer case at half-filling. Of course, there is 
no constraint that hn = hi%, but we suspect that they 
are nevertheless sufficiently correlated to reduce the sign 
problem. 

V. CONCLUSIONS 

In this paper we have used DQMC simulations to de- 
termine the phase diagram, in the (tj_/t, U) plane, of the 
half-filled bilaycr Hubbard model. Our phase diagram 
exhibits metallic, band insulating and Mott insulating 
phases in qualitative agreement with CDMFT rcsultsAS,. 
However, the entire metallic phase we find is paramag- 
netic with no antiferromagnetic metallic regions. 

In addition, we have shown that the doped bilayer 
Hubbard Hamiltonian has an attractive d-wave pairing 
vertex which approaches close to TdPd = —1, signaling 
a superconducting transition. This value is much more 
singular than that observed in the single layer model, 
partly because it is more attractive when compared at 
the same inverse temperature, and partly because it is 
possible to simulate to values of (3 which are two to three 
times larger than for a single plane. However, the un- 
corrected Pd gets small, so that the enhancement of the 



correlated Pd is not very dramatic. On the other hand, 
and unlike what happens in the single layer d = 2 Hub- 
bard model, the enhancement here is maximum when the 
system is doped, in agreement with the phenomenology 
of cuprate superconductors. 

Pairing in systems with separate electron and hole 
doped sheets has a long history in the context of ex- 
citon condensation^, but our primary motivation here 
has been the recent report of cuprate-based systems 
Ba2Ca3Cu40gF2 and Ba2Ca2Cu40eF2 which have both 
types of doping*^. Our results for the magnetic and 
pairing correlations bear interesting connections to those 
materials. While the bilayer simulations reported here 
contain the the essential feature of coupled electron and 
hole doped layers, it is natural to consider direct numer- 
ics of three and four layer compounds. Such studies will 
require an order of magnitude greater simulation time, 
and also have an at present unknown sign problem. In 
general the sign problem becomes worse with lattice size 
(and hence with number of layers), but it is tempting to 
speculate that, if the above picture of correlated determi- 
nant signs presented above is correct, similar correlations 
might provide protection for larger numbers of layers. 
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